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A method for obtaining a statistically 
stationary turbulent free shear flow 

By S. F. Timson, S. K. Lele AND R. D. Moser 
1. Motivation and objectives 

The long-term goal of the current research is the study of Large-Eddy Simulation 
(LES) as a tool for aeroacoustics. New algorithms and developments in computer 
hardware are making possible a new generation of tools for aeroacoustic predictions, 
which rely on the physics of the flow rather than empirical knowledge. LES, in con- 
junction with an acoustic analogy (Lighthill 1952), holds the promise of predicting 
the statistics of noise radiated to the far-field of a turbulent flow. 

While there have been preliminary studies where LES was used in aeroacoustic 
calculations, a thorough examination of LES’s predictive capabilities has not been 
undertaken. Recent advances in subgrid-scale models (Ghosal ei al. 1992) have 
shown promising results, but accurate acoustic predictions present a far more strin- 
gent test of LES’s fidelity. LES does not resolve all scales in the simulated flow, and 
results of a given simulation depend on the subgrid-scale model employed. Previous 
verification of LES results has focused on one point statistics. The application of 
an acoustic analogy requires a two point space-time correlation, and it is unclear 
how accurately LES will reproduce this quantity. 

It has been shown that the dominant features of far-field noise are associated with 
the energy containing range of scales in the near-field, offering hope that LES can 
supply accurate predictions. However, in practical applications, the subgrid-scale 
energy can represent a significant fraction of the total. Further, as Crighton (1988) 
points out, care must be taken when applying acoustic analogies. Computation- 
ally, spurious but efficient low order acoustic sources may result from discretization 
errors, boundary conditions, etc. It is important to assess whether or not the 
subgrid-scale model acts as such a low order acoustic source. 

LES’s predictive ability will be tested through extensive comparison of acoustic 
predictions based on a Direct Numerical Simulation (DNS) and LES of the same 
flow, as well as a priori testing of DNS results. The method presented here is aimed 
at allowing simulation of a turbulent flow field that is both simple and amenable to 
acoustic predictions. A free shear flow that is homogeneous in both the stream wise 
and spanwise directions and which is statistically stationary will be simulated using 
equations based on the Navier Stokes equations with a small number of added terms. 
Studying a free shear flow eliminates the need to consider flow-surface interactions 
as an acoustic source. The homogeneous directions and the flow’s statistically sta- 
tionary nature greatly simplify the application of an acoustic analogy. 

2* Accomplishments 

A method allowing simulation of a statistically stationary free shear flow has been 
developed. The method is an extension of that presented by Spalart (1988) to the 
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case of a wake or coflowing jet in the small deficit limit. The results are similar to 
those for the sink-flow boundary-layer presented in Spalart (1986). The derivation is 
more rigorous than either of the above analyses due to the simplification introduced 
by explicitly considering the small deficit limit. Some limited testing of the method 
has been carried out, and a more detailed validation is in progress. 


2.1 Mathematical formulation 

The flow to be simulated is that of a plane wake or coflowing jet. Self similar 
behavior is assumed and forms the basis for the rest of the analysis. The classical 
conditions necessary for self-similarity in a plane wake (in the small deficit limit) 
are 
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where 6 is the wake thickness, U 0 = Au max is the maximum velocity deficit, and x 
denotes the streamwise direction. 

A coordinate transformation is defined such that 


(x,y,z,t) -i (x,r],z,t) 


( 3 ) 


where tj = y/8(x). The spatial coordinates have been normalized by some initial 
length scale L 0 , and time is normalized by some initial time scale t a . Thus all 
quantities are non-dimensional. In the new set of coordinates, the profile thickness 
is independent of x. Lines of constant tj have slope S where for the above coordinate 
system 

e d8 

s ="j- x - < 4 > 

The Jacobian of the coordinate transformation in space may then be written 
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where T = 6 for the coordinate change given in (3). The actual Jacobian is a 
4x4 matrix, but as no transformation is carried out in time, the extra elements are 
deleted for simplicity. 

Following Spalart (1988), a transformation is made to the dependent variables as 
well. Let the Cartesian velocity components be denoted by (u*, v*, w*), then the 
contravariant velocity components will be denoted by (u, v, w). The two sets are 
related through J ~ l , giving 
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Use of the contravariant velocity components preserves the form of the transport 
terms in the new coordinate system, i.e. 


*d m d m d _ d J , J 

u +v + w = u x~ + 

ox oy uz ox or} oz 

As Spalart shows, the continuity equation becomes 
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The x-momentum equation becomes 
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where p is the kinematic pressure. 

Normalized velocity components are introduced, and the equations are expressed 
in terms of the velocity deficit. The variable change is 


/ u(z,r/,z,f) \ 

v(x,T},Z,t) \ 
w{x,r],z,t) I 
\p(x,rj,z,t) ) 
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( 10 ) 


The maximum deficit is used to scale the velocity components so that the mean 
and r.m.s of the velocity components are independent of downstream location in 
the new coordinate system. 

The resulting equations are not presented as they are quite cumbersome. One 
should note that new terms involving x derivatives of the velocity scale are intro- 
duced. These terms are recast in terms of the undifferentiated velocity deficit using 
the assumed self-similar streamwise evolution. Using Eqs. (2) and (4) and recalling 
that T = <5, it can be shown that 
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and 



The equations are further simplified by choosing the length scale L a to be <5 0 , the 
dimensional value of the layer thickness at the downstream station to be simulated. 
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This implies the dimensionless quantity 6 is 1. Then T = 1, rj derivatives of T are 
zero, and x derivatives are removed through the identity 


F x — Sij. 


Further, if 6 is unity, rj and y become equivalent, and all instances of rj in the 
equations may be replaced with y. Then defining the parameter 


the continuity equation becomes 


t* U X + Vy + W z = 0. 

€ * 


The x-momentum equation becomes 

il t + UUx + Wy + WU Z 

( U X *5* y U 2^ - 

— SyU J = —p x -f Sp y + viscous terms. (16) 

The viscous terms are quite complicated and are omitted for simplicity. 

A Galilean transformation is now applied. A new coordinate x ' is defined such 
that 

a;' = x-Uoot. (17) 

Applying this transformation, the terms in the momentum equations resulting from 
convection at U^o (e.g. the u x /e term in (16)) vanish. 

The slope of the new coordinate lines 5 will now be expressed in the new coor- 
dinate system. Returning to the unmoving coordinate system, using (2) and (4) it 
can be shown that 

s- v 44-AhV m 
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A quantity that was only a function of x in the stationary frame becomes a function 
of both x ' and t in the moving frame, i.e. 

F(x) = F(x f + U OQ t). (19) 

Differentiating with respect to time in the new frame shows 

1 dF dF 

Uoo dt ~ ~ dx' 

Applying (20) to (18) gives S in the new coordinate system as 

s = ¥1 Q?) e » ' < 21 > 
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A parameter a is defined by 


U„S d ( 1 \ 

2 dt\U l)' 


( 22 ) 


and the scaling of U 0 in (2) implies that a is constant. It follows that S = aey. 

The small deficit limit is then enforced, letting e->0, giving for the continuity 
and momentum equations 

V •« = -<*, (23) 


and 


u t + uj x u — — Vp + a 
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(24) 


These are the equations that will be integrated to yield the turbulent flow field. Note 
that (20) indicates that V = - ay , and that the fluctuating velocity components 
are divergence free. Also note that all added viscous terms dropped out when the 
limit was taken. 


2.2 Implementation 

The modified equation set was implemented in the spectral code used by Rogers 
and Moser (1993) in their study of self-similar turbulent mixing layers. The nu- 
merical method (Spalart, Moser and Rogers 1991) employs Fourier analysis in con- 
junction with periodic boundary conditions in both the streamwise and spanwise 
directions. The application of periodic conditions in the streamwise direction is an 
approximation as the length and time scales for the fluctuating quantities are x de- 
pendent. The approximation should be a good one for the normalized components 
in this coordinate system as the variation of these scales is slow compared with the 
variation of 6. The non-homogeneous direction is handled through the use of rapidly 
decaying spectral basis functions in conjunction with two slowly decaying additional 
functions that exactly represent the irrotational component of the solution far from 


the vortical region. 

The major question in implementing the modified equations is how to evaluate a. 
An initial value may be calculated for a by taking the second moment of the mean 
x-momentum equation, setting Ut to zero, and solving for a. Making use of the fact 
that the mean field is a function only of y, and integrating by parts, discarding the 
boundary terms, gives 


/ ^y 2 (J x J) r + j^Udy 

Jy 2 Udy 


(25) 


This value is then held, and the mean flow is allowed to evolve until a statistical 
steady state is reached. The flow in the interim has no physical meaning. 

It may appear that a similar scheme may be used to constantly update a as 
the flow is evolving, in effect setting a to the value it would have if the flow was 
stationary at that moment, hopefully speeding evolution to the stationary state. In 
practice, however, this presents problems. It appears that a feedback loop involving 
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the Reynolds stresses results. Further investigation in this area is clearly needed. 
For the present work this is not critical as the initial conditions to be used are fully 
developed wakes simulated by Moser and Rogers (1994). Calculating a based on the 
initial condition and fixing it should not result in a significant evolutionary period. 

2.3 Testing 

The first test case was a laminar wake. The initial condition was simply a Gaus- 
sian. The flow evolved into a shape that was nearly Gaussian (a Gaussian does not 
solve the modified equation set) and remained stationary. A second test consisted 
of an initial condition of y 2 times a Gaussian, such that the profile had the same 
mass flow and second moment. The flow evolved to the same steady solution as was 
reached starting with a Gaussian. 

The first turbulent test case is that of a low Reynolds number turbulent wake. 
This test is currently underway. The initial condition was taken from a temporally 
evolving simulation done by Moser and Rogers (1994). The results will be compared 
qualitatively to the results obtained from the temporal simulation to ensure that 
the results produced by the modified equation set are reasonable. 

3. Future plans 

Upon completion of the testing with the low Reynolds number wake, a higher 
Reynolds number DNS will be carried out, building the database necessary for 
aeroacoustic predictions. A priori testing of the DNS results will be carried out, 
comparing the acoustic prediction from the DNS data and a filtered version of the 
same flow-field. It has been seen previously that LES simulations Eire capable of 
producing better results than indicated by a priori tests. Therefore, the code will 
also be modified, and an LES of the same flow will be carried out. Using these 
databases, the issues discussed in section 1 will be investigated. 
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